Thermal Rectification in Graded Materials 
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In order to identify the basic conditions for thermal rectification we investigate a simple model 
with non-uniform, graded mass distribution. The existence of thermal rectification is theoretically 
predicted and numerically confirmed, suggesting that thermal rectification is a typical occurrence 
in graded systems, which are likely to be natural candidates for the actual fabrication of thermal 
diodes. In view of practical implications, the dependence of rectification on the asymmetry and 
system's size is studied. 
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The study of the underlying dynamical mechanisms 
which determine the macroscopic properties of heat con- 
duction has opened the fascinating possibility to control 
the heat current. In particular a model of a thermal 
rectifier has been proposed [l[ and since then, the phe- 
nomenon of thermal rectification has been intensively in- 
vestigated 0-Q in order to analyze and improve the rec- 
tification effect, including experimental realizations (Io| . 

However, as correctly pointed out in Ref. [3[, most 
recurrent proposals of a thermal diode, based on the se- 
quential coupling of two or three segments with different 
anharmonic potentials, are difficult to be experimentally 
implemented and the rectification power typically decays 
to zero with increasing the system size. In addition, most 
investigations so far have been based on numerical sim- 
ulations and a much better theoretical understanding is 
highly desirable both for fundamental reasons as well as 
for obtaining useful hints for the actual realization of de- 
vices with satisfactory rectification power. 

Along these lines, graded materials are attracting more 
and more interest: Papers with numerical 0^3; analyt- 
ical 0, Q and even experimental [Io| studies have ap- 
peared recently in the literature. 

The present paper addresses the fundamental dynam- 
ical mechanisms that lead to rectification. Our strategy 
is to consider a simple model that contains the minimal 
ingredients we theoretically judge to be necessary to rec- 
tification, and compare the numerical results with the 
theoretical predictions. As the features of our model are 
also shared by more realistic models such as anharmonic 
chains of oscillators, we conjecture that the obtained re- 
sults may have practical implications as well. Our study 
allows us to understand the basic ingredients behind rec- 
tification, to describe non-trivial and important proper- 
ties of the heat flow, and to show that rectification in 
graded materials could be a ubiquitous phenomenon. 

We consider a chain ll| of elastically colliding parti- 



FIG. 1: (Color online) The schematic plot of our model. Dot- 
ted lines divide elementary unit cells. In each cell there is 
a bar which is subjected to elastic collisions with both cell 
boundaries and with neighboring bullets. The first (last) cell 
is coupled to a heat bath at temperature tl (tr). 



and "bullets" , respectively. (See Fig. 1.) Each bar is con- 
fined inside a cell of unit length; that is, besides elastic 
collisions with its neighboring bullet(s), it is also subject 
to elastic collisions with the two boundaries of the cell. 
A bullet only undergoes collisions with its two neigh- 
boring bars. Apart from collisions, all particles move 
freely. We denote the masses (velocities) of bars by Mi 
(vi), I — 1, ••• ,Z, and those of bullets by nik (wfc), 
k = 1, • • • ,Z — 1, respectively, with Z being the total 
number of cells. Two statistical thermal baths with dif- 
ferent temperatures tx and tr are put into contact with 
the two ends of the system: When the first (last) bar 
collides with the left (right) side of the first (last) cell, it 
is injected back with a new velocity determined by the 
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cles of two kinds referred to, in the following, as "bars" 



The Boltzmann constant ks is set to be unity. In our 
molecular dynamics simulations, after the system reaches 
the steady state, we can compute the local temperature 
of each cell and the steady heat flux across the system. 
The local temperature of the Zth cell is defined as the ki- 
netic temperature of the bar, i.e., r; = (Mivf /ks), where 
(•) stands for the time average. The steady heat flux is 
measured as the time average of the energy exchanged 
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in unit time between the left heat bath (the last bar) 
and the first bar (the right heat bath), or that between 
any two neighboring particles. We have verified that the 
so-measured heat flux is the same as the local heat flux 
carried by each particle, which is defined as {Mivf /2) 
((miuf/2)) for the Ith bar (bullet). 

In order to investigate rectification, we consider the 
mass graded version of the above chain. Precisely, we 
will focus on a system with mass graded bullets m\ < 
TO2 < • • • < mz—i and identical bars, Mi = ■ ■ ■ = Mz = 
1. In order to avoid boundary effects due to coupling 
with the heat baths, the rectification effects are measured 
numerically only over the N central cells, with the first 
(last) Zl {Zr) not being considered. Hence the effective 
system size is N = Z — Zl — Zr. 

The starting point of our theoretical approach is the 
expression for the local heat flow inferred from the ho- 
mogeneous case. We recall that the Fourier law of ther- 
mal conduction has been shown numerically to hold in 
the homogeneous version (i.e., the bullets also have the 



same mass m; = m) of our model We also recall 

that homogeneous lattice models of oscillators with, e.g., 
quartic anharmonic on-site potentials, have been stud- 
ied by different analytical methods 13 - 1(| and shown to 
follow the Fourier law as well. It implies that for these 
systems, in the steady state the local heat current at 
position x is proportional to the temperature gradient 
F(x) — —JC(x)VT(x), where the heat conductivity K(x) 
is a function of the local temperature T(x) (and of other 
parameters). Precisely, for the homogeneous models of 
anharmonic oscillators, such a local expression reads 



F=(T j -T j+1 )/Ct 
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where Tj is the local temperature at the jth site and 
Tj 3 = (T/+T/ +1 )/2 with p w 2. For the theoretical study 
of rectification in our inhomogeneous chain, we start from 
the inhomogeneous version of such equation, i.e., we write 
the heat flow from cell j to j + 1 as 



Tjj+i = K,(VI")j 
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Note that the homogeneous version of our model obeys 
the Fourier law with a thermal conductivity that scales as 
T 1 / 2 (see the following), in contrast to the 1/T 2 behavior 
of the above anharmonic oscillators models. 
For the heat flow in the steady state we have 
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Hence, by summing up the (N — 1) equations in ([3]) for 
j = i,2,-.. ,N- 1, we get T = K{T X - T N )/(N - 1), 
where 
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FIG. 2: (Color online) The numerically computed heat con- 
ductivity for the homogeneous chain of size ./V = 45 at tem- 
perature T = 1 (see Ref. [HJ). The bars have unit mass and 
the bullets have the same mass m. The two solid (dotted) 
lines indicate the bullet mass region used for the simulations 
in Fig. 3 (b) [Fig. 3(c)]. 



Then the Fourier law follows in the case where the ther- 
mal conductivity K, remains finite when N —¥ oo. 
From Eqs. © and (gj), it follows that 

(Ti - TaXdTf + C 2 T 2 Q ) = (T 2 - T 3 )(C 2 T 2 Q + C 3 T 3 Q ) 
= • • • = (Tn-i — Tn)(Cn-\Tn-i + CnT^); (6) 

Thus, given the temperatures of the first and the last cell, 
i.e., Ti and T/v, by using the equations above we may de- 
termine the inner temperatures T2, T3, • • • , Tjy-i. As it 
may be a hard problem to obtain the analytical solution 
of these equations, we may turn to numerical computa- 
tions, or else, as in Refs. [1,0], we may assume small tem- 
perature gradients in order to simplify the computations. 
In this latter regime, i.e., for Ti = T+a\e, Tjy = T+a^e, 
e small, we have T& = T + a^e + C(e 2 ), where have 
to be determined [we carry out the computations only 
up to 0(e)]. Algebraic manipulations give us a% = ai + 
{a N - ai )C{k)/C{N), for k = 2,--- ,N-1, where C{k) = 
[(Ci + C2)- 1 + (C 2 + Cs)- 1 + ■ ■ • + (C fc _i + Cfe)- 1 ] . The 
thermal conductivity defined in Eq. ([5]) is therefore 
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Now we analyze the heat flow for the system with in- 
verted thermal baths. In this case we denote the tem- 
perature of the jth cell by Tj and assume that T[ — 



Tjy and T' N — 
k = 2,3, ••■ ,N 
a' k = a N - (a N - 



Ti. 



We also write T' k — 



T + ale for 



- 1, and after the computations we get 
ai)C(k)/C(N). Recalling that a[ = a N 
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FIG. 3: (Color online) (a) Temperature 
profile of the graded chain with bullets' 
masses increasing in the interval (1,2) for 
tl = 1.2 and tr — 0.8 (red dots), and 
inversely, t' l — 0.8 and t' r — 1.2 (blue 
squares). The dashed lines indicate the 
end temperatures Ti (Tn) and T N (Tl) 
of the central segment considered in our 
computations, (b)-(d) Heat conductiv- 
ity measured over the central segment 
of the graded chain with (b) and (c) for 
bullets' masses increasing in the interval 
(1,2) and (0.4,0.8), respectively (all bars 
have unit mass), and (d) for the graded 
chain where all the masses of bars and 
bullets increase in the interval (1,2). See 
text for details. 



and a' N = a±, one obtains for the thermal conductivity 
K! of the "inverted" system, an expression similar to JC 
(with a'j replacing Oj). We then have: 
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The homogeneous model, where all bars have a unit 



mass and all bullets have mass m, has been shown [11 1 
to follow the Fourier law and the thermal conductiv- 
ity ft(m, T) at temperature T = 1 is as shown in Fig. 
2. As the dynamics of our model can be essentially 
described by a series of particle collisions, which is in- 
variant under a time rescaling t — > t/7 (with 7 being a 
positive constant), or equivalently a temperature rescal- 
ing T — ¥ y^T, we have «;(m, T) = K,(m, l)VT. There- 
fore in the limit of small temperature gradients, Eq. ^ 
leads to a = 1/2. To conclude our theoretical analy- 
sis we still have to determine the coefficients Cj. The 
behavior of these coefficients as function of m can be de- 
duced from Fig. 2 for the homogeneous chain. Indeed 
let us also assume for our inhomogeneous system a small 
mass gradient and determine such coefficients by prop- 
erly relating them to their homogeneous versions; i.e., 

Cj = [K(m j+ZL -i, 1) + K(m,j + z L ,l)]/2- 

The qualitative picture which now emerges is quite 
clear: From Eq. (jSJ), in particular from the difference be- 
tween the local thermal conductivities, we predict, and 
will later numerically confirm, the properties of rectifi- 
cation. First, from Eq. (|SJ) it immediately follows the 
existence of rectification which is particularly clear when 
/t(mj,l), and so Cj has a monotonic behavior. It also 
follows that rectification increases with the gradient of 
temperature. Moreover, from Eq. ([5]). one can predict a 



larger flow from a smaller to larger mass direction in the 
region where K(rrij, 1) decreases with rrij, and the oppo- 
site behavior in the region where n(m,j, 1) increases with 
rrij. The fact that our approach allows to predict not 
only rectification of the heat current but also its direc- 
tion is quite interesting. In this respect we recall that 
rectification in a graded system with the heat current 
larger in the direction of decreasing particle masses has 
been observed experimentally in nanotubes with nonho- 
mogeneous external mass loading [10] and numerically in 
the mass graded Fermi-Pasta-Ulam model 0. Instead, 
a case with the heat current larger in the reversed direc- 
tion (of increasing masses) has been found by molecular 
dynamics simulations in mass graded ideal gases @ and 
some carbon nanotubes Q. 

We now turn to the molecular dynamics studies of our 
system which not only confirm our theoretical predictions 
but also extend, in fact, our analytical results to regimes 
beyond small gradients of temperature and mass, where 
the rectification phenomena become more relevant. To 
this end let us consider first a chain of Z = 32 cells 
with the bullets' masses increasing in the range (1,2) 
(i.e. the interval between the two solid vertical lines in 
Fig. 2) with mi = 1 + l/Z . After the steady state 
is reached we calculate the time averaged temperature 
profile and the heat current. Averages are taken over 
a number of collisions per particle larger than 2 x 10 9 . 
In Fig. 3(a) we show the computed temperature profile 
for tl = 1.2 and tr = 0.8, as well as for t' l = 0.8 and 
t' r — 1.2. To evaluate the heat conductivity K, and K' 
we neglect the first (last) Zl {Zr) cells in order to get 
rid of the boundary effects. The values of Zl and Zr 
are determined by ensuring that T\ — T' N and Tn = T[. 
(Note that T x = t 1+Zl and T[ = t' 1+Zl for I = 1, • ■ ■ , N.) 
In this case we have Zl = 5 and Zr = 2. [See Fig. 
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3(a).] The value of thermal conductivity JC — F(N — 
l)/(Ti — Tjv) is then calculated over the central segment 
of N = 25 cells with numerically measured J 7 , Ti and 
TV- KJ is computed in the same way. 

The dependence of the heat conductivity JC and JC' on 
AT for T = 1 is given in Fig. 3(b). Here T\ — T' N — 
T + AT and T N = T[ = T - AT. It is clearly seen that 
a system with graded mass particles does present rectifi- 
cation, and that rectification increases with the temper- 
ature gradient. In addition, as in this case Cj decreases 
monotonically (see the region between two solid vertical 
lines in Fig. 2), according to our theoretical analysis the 
value of JC for T\ > Tv should be larger than JC' for 
T[ < T' N . This is confirmed by the numerical results. 

One advantage of our model is that, as shown in Fig. 2, 
it can display regions where Cj decreases with j as well as 
regions where Cj increases with j as, e.g., in the interval 
between the two dashed vertical lines in Fig. 2. In the 
latter situation we predict that the behavior of JC and JC' 
will be opposite to the previous case, i.e., that JC for T\ > 
Tjv will be smaller than JC' for T[ < T' N . This is confirmed 
by Fig. 3(c) where we present the numerical results for 
a graded chain with the masses of bullets in the interval 
(0.4,0,8) with mi = 0.4(1 + l/Z). All the computations 
are done as before with Z = 32, Zl = 3, Zr = 4, N = 25 
and T = 1. 

An obvious interesting question is how to enhance rec- 
tification in our graded chain. As our theoretical analysis 
suggests, a possible way to increase rectification is to in- 
crease the asymmetry of the system. We have therefore 
analyzed the case where the masses of the bars are also 
allowed to have a graded distribution. To be precise, 
all the masses are graded in between 1 and 2. Namely 
Mi = 1 + (Z-0.5)/Zfor I = l,-- - , Zand to/ = 1 + Z/Zfor 
I = !,-■■ , Z - 1. Here we take Z = 32, Z L = 5, Z R = 1, 
N = 26 and T = 1. The results are depicted in Fig. 
3(d) and show a significant increase in rectification. We 
stress that the increase of rectification with asymmetry 
is an important effect which may be of practical inter- 
est: Similar mechanisms may be present in more realistic 
models due to graded interparticle interactions and/or 
graded on-site potentials, in addition to graded masses. 

Finally, we remark that for the case presented in Fig. 
3(d) the heat conductivity K. and JC' essentially do not 
change with the system size. This has been verified by 
additional simulations with Z = 64, 128 and 256 (with 
Z L = 8, 15, and 29, Z R = 2, 4, and 8, and N = 54, 109, 
and 229, respectively). The rectification, defined as 

= (JC — JC')/ JC' , does not change with the system size 
as clearly illustrated in Fig. |4] The fact that rectifica- 
tion does not vanish with increasing the system size is an 
important property, which is absent in the known diode 
models given by the sequential coupling of different seg- 
ments. 

In summary, in this work, we have performed analyt- 
ical and numerical investigations of the heat flow in a 
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FIG. 4: (Color online) The rectification measured over the 
central segment of the graded chain where all the masses of 
bars and bullets increase linearly in the interval (1,2). The 
dots, squares, triangles and stars represent the results for sys- 
tem size Z = 32, 64, 128 and 256, respectively. 



simple model, devoid of intricate interactions, in order to 
find results of general validity. Our results demonstrate 
that thermal rectification takes place in asymmetric sys- 
tems with local heat flow proportional to the tempera- 
ture gradient and with the local conductivity depending 
on temperature. This conclusion is also supported by ad- 
ditional numerical investigations 17], along the lines of 



Refs. [Ill and [18j , of the energy diffusive behavior and 
of correlation properties of energy fluctuations. Finally, 
it is interesting to notice that in our model the diffusive 
behavior (typical of the Fourier law) is associated with 
the thermal rectification phenomenon. 
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